package scu.maqiang.bes;

import java.util.Arrays;

import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.T3Type;
import scu.maqiang.numeric.DirectDSolver;
import scu.maqiang.numeric.MVO;

public class PossionBE2DT3L0 {

	public static void main(String[] args) {
		Mesh2T3 mesh = new Mesh2T3().square2D(20, 20, T3Type.Right);
		//mesh.displayBoundarys();
		BES2DT3L20 bs = new BES2DT3L20(mesh);
		//bs.displayMiddlePoint();
		bs.setDirichletBoundaryValue(0, 1, 2, 3, 4);
		double[][] H = new double[bs.GetNdof()][bs.GetNdof()];
		double[][] G = new double[bs.GetNdof()][bs.GetNdof()];
		bs.assembleLaplaceMatrix(new double[] {1.0}, H, G);
		double[] source = new double[bs.GetNdof()];
		bs.assembleSource(new double[] {1.0}, source);
		//System.out.println(MVO.toString(source));
		double[] knowValue = bs.getKnownValue();
		//System.out.println(MVO.toString(knowValue));
		double[] RHS = MVO.matmul(G, knowValue);
		//System.out.println(MVO.toString(RHS));
		Arrays.setAll(RHS, i -> RHS[i] + source[i]);
		
		//System.out.println(MVO.toString(knowValue));
		DirectDSolver solver = new DirectDSolver(H);
		double[] x = new double[bs.GetNdof()];
		solver.PLU();
		solver.PLUSolve(RHS, x);
//		System.out.println(Arrays.toString(x));
		double[][] result = bs.getSolution(x);
		//System.out.println(MVO.toString(result));
		double[] constSource = {1.0};
		//double[][] result2 = bs.computeDomainValues(result[0], result[1]);
		System.out.println(bs.computeValue(result[0], result[1], constSource, new double[] {0.25, 0.25}));
		System.out.println(bs.computeValue(result[0], result[1], constSource, new double[] {0.50, 0.25}));
		System.out.println(bs.computeValue(result[0], result[1], constSource, new double[] {0.75, 0.25}));
		System.out.println(bs.computeValue(result[0], result[1], constSource, new double[] {0.25, 0.50}));
		System.out.println(bs.computeValue(result[0], result[1], constSource, new double[] {0.50, 0.50}));
		System.out.println(bs.computeValue(result[0], result[1], constSource, new double[] {0.75, 0.50}));
		System.out.println(bs.computeValue(result[0], result[1], constSource, new double[] {0.25, 0.75}));
		System.out.println(bs.computeValue(result[0], result[1], constSource, new double[] {0.50, 0.75}));
		System.out.println(bs.computeValue(result[0], result[1], constSource, new double[] {0.75, 0.75}));
		//System.out.println(bs.computeValue(result[0], result[1], constSource, new double[] {1, 1}));
		
		double[][] domainResult = bs.computeDomainValues(result[0], result[1], constSource);
		mesh.toTecplot("BEMPossionResult.dat", domainResult);
	}

}
